Particle-hole symmetry in a sandpile model 
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In a sandpile model addition of a hole is denned as the removal of a grain from the sandpile. We 
show that hole avalanches can be defined very similar to particle avalanches. A combined particle- 
hole sandpile model is then defined where particle avalanches are created with probability p and 
hole avalanches are created with the probability 1 — p. It is observed that the system is critical 
with respect to either particle or hole avalanches for all values of p except at the symmetric point of 
p c = 1/2. However at p c the fluctuating mass density is having non-trivial correlations characterized 
by 1// type of power spectrum. 
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The dynamics of a large number of physical processes 
are characterized by bursts of activity in the form of 
avalanches. For example, the mechanical energy release 
during earthquakes Q , river networks Q , forest fires Q , 
land slides on mountains or sand avalanches on sandpiles 
etc. Bak, Tang and Wiesenfcld in 1987 proposed that 
these systems may actually be exhibiting signatures of 
a critical stationary state More precisely they sug- 
gested that as an indication of the critical state, long 
ranged spatio-temporal correlations may emerge in some 
systems governed by a self-organizing dynamics, in ab- 
sence of a fine tuning parameter. This is in essence the 
basic idea of Self-organized Criticality (SOC). Sandpile 
models are the prototype models of SOC |1 IE II 0, Ef El • 

In the sandpile model an integer variable n represent- 
ing the number of particles (sand grains) in a sand col- 
umn is associated with every site of a square lattice. The 
system is driven by adding a particle at a randomly se- 
lected site i: rii — > m + 1 . A threshold number n v c for the 
stability of a sand column is pre-assigned. If at any site 
the particle number n% > nP c the column topples and this 
site looses 4 particles and all four neighboring sites j get 
one particle each p|. 

TT-i — s- — 4 and rij — ► rij + 1 (1) 

As a result some of these neighboring sites may also top- 
ple which creates an avalanche of sand column topplings. 
The extent of such cascading activity measures the size 
of the avalanche. 

Sand particles drop out of the system through the 
boundary of the lattice so that in the steady state the 
fluxes of in-flowing and out-flowing particle currents bal- 
ance. In a stable state number of particles at all sites are 
less than n§. Addition of a particle takes the system from 
one stable state to another stable state. Dhar had shown 
that under this sandpile dynamics, a system evolves to 
a stationary state where all stable states are restricted 
to a subset of all possible stable states. These states are 
called recurrent states and they are characterize d by the 
absence of forbidden sub-configurations (FSCs) |lflj |. All 
recurrent states occur with uniform probabilities in the 
stationary state. A stable state which is not a recurrent 
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FIG. 1: The fluctuation of the mean number of particles per 
site with time in a system of size L — 64 for the probabili- 
ties p = 0.60 (top) 0.52 (middle) and 0.50 (bottom). Both 
the width of fluctuation and the correlation increases as p 
approaches p c — 1/2. 



state is called a transient state and never appears in the 
stationary state. 

The main question we would like to ask in this paper 
is, for an arbitrary sandpile model to attain the BTW 
critical behaviour is it absolutely necessary that the sta- 
tionary states should only be the recurrent states of the 
BTW model? Can it happen that the neighbouring tran- 
sient states which are very close to the recurrent states 
of the BTW model are also acceptable in the stationary 
states to achieve the BTW critical behaviour? In the fol- 
lowing we introduce the concept of holes and on addition 
of holes to the system the stationary states of the result- 
ing sandpile model cannot be anymore strictly restricted 
to the recurrent states of the BTW model since the FSCs 
can very well be present in the stationary states of this 
model. In fact the recurrent and stationary sets coin- 
cide in this model. However the distribution of weights 
of these states may be quite non-trivial and this question 
remains open. Our numerical results show that even in 
such a case the critical behaviour is very similar to that 
of the BTW model. 
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FIG. 2: The time averaged number of particles per site 
(n(p,L)) in the stable stationary states as a function of the 
probability of adding a particle p. This variation is symmet- 
ric about the mid point p c = 1/2 and (n(p c ,L)) = 3/2. The 
data is for L = 32 (circle), L = 64 (square) and for L = 128 
(triangle). 

A 'hole' may be defined as the absence of a particle. 
Therefore adding a hole to a lattice site implies taking 
out one particle from that site: rii — » — 1. Repeated 
addition of holes at randomly selected sites may reduce 
the number of particles at a site less than another pre- 
assigned threshold n^. Therefore if rii at a site is less 
than n c l , the site losses four holes i.e., four particles are 
added to this site and each neighboring sites gets one 
hole (looses one particle): 

m — > rii + 4 and rij — > rij — 1 (2) 

We call this event as a 'reverse toppling'. Consequently 
at some of the neighboring sites particle numbers may 
also go below the nj? which again reverse topple and thus 
an avalanche of reverse topplings take place in the sys- 
tem. Addition of a particle creates a particle avalanche 
where as the addition of a hole creates a hole avalanche. 
We assign n v c = 3 and n£ = 0. 

Inverse avalanches were introduced before to get back 
the recursive configuration corresponding to the particle 
deletion operator [13j. 

During the particle avalanches particle current flows 
into the system by addition of particles in the bulk of the 
system and then they flow out of the system through the 
boundary. On the contrary in hole avalanches particle 
current flows into the system through the boundary and 
flows out of the system through the bulk of the system. 
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FIG. 3: Variation of the correlation time r(p,L) of fluctu- 
ation of the particle density (n(p,L)} with the probability p 
for system size L. The data is shown for two system sizes: 
L = 32 and 64. 



In this paper we study a combination of particle and 
hole avalanches. We probabilistically add either a parti- 
cle with a probability p or add a hole with a probability 
1 — p. Therefore when p = 1, the situation is identical 
to the ordinary BTW model of sand avalanches when no 
hole is added. On the other hand for p = only holes 
are added to the system and no particle. Therefore for 
p > 1/2 more particles are added to the system than 
the number of holes and therefore the net particle cur- 
rent is 2p — 1 and it flows from the bulk of the system 
to the boundary. However for p < 1/2 holes are dropped 
more than particles and the net particle current flows in 
the opposite direction. At p = 1/2 there is no net cur- 
rent in the system. We consider p = p c = 1/2 is a critical 
probability and study the behavior of this system around 
this critical probability. The time t is measured by the 
number of particles and holes dropped in the system. 

If rii(t,p,L) is the generalised notation for the num- 
ber of particles at site i, then the total number of par- 
ticles in the system is: n(t,p,L) = T,f =1 rii(t,p, L). The 
mean number of particles per site is then (n(t,p,L)) = 
n(t,p, L)/L 2 . In the stationary state (n(t,p,L)) fluctu- 
ates rapidly around its time averaged value (n(p, L)). 
These fluctuations are shown in Fig. 1 for the system 
size L = 64 and for p = 0.60, 0.52 and for 0.50. It is 
observed from this figure that both the width as well as 
the correlation of fluctuation increases as p approaches 
p c from either side of it. 

The time averaged number of particles per site 
{n(p, L)) is a function of p and the system size L. At 
p = 1 it is equal to the average number of particles per 
site in the ordinary BTW model which is ni = 2.125 in 
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FIG. 4: (a) The autocorrelation C(t,p c , L) of the time series of 
the fluctuating mean number of particles per site {n(t,p c ,L)) 
for the system size L = 64 at p c — 1/2. C(t,p c ,L) is plotted 
with the scaled time axis t/I 2 on a semi- log scale. In (b) the 
power spectrum S(f,p c ,L) is plotted with the frequency / 
on a double logarithmic scale showing a power law decay of 
the spectrum with the spectra exponent being nearly equal 
to one. 



the asymptotic limit of large system sizes jjj ^| . As p 
decreases (n(p,L)) slowly decreases but near p c = 1/2 it 
decreases very fast to a value of (n(p, L)) —3/2. When p 
decreases from 1/2 even further, (n(p,L)) decreases fast 
but eventually saturates to a value of uq = 3— n\ = 0.875. 
In Fig. 2(a) we show this variation. To see if the steep 
rise of {n(p, L)) around p c is associated with some critical 
exponent, we make a scaling plot of (n{p, L)) with p — p c 
for a number of different system sizes L in Fig. 2(b). The 
data collapse shows: 



<n(p,L)>~S((p-p c )L 



1.68\ 



(3) 



The width of fluctuation is calculated as: w(p, L) = 
(n 2 (p,L)) — (n(p,L)) 2 . For a given L the width is max- 
imum at p = p c and then monotonically decreases as 
\p — p c \ increases. On the other hand for a given p the 
width also decreases with increasing L. It is observed 
from numerical estimation that at p = p c , w(p c , L) de- 
creases with system size as w(p c ,L) — w a + wiL^ 1 / 2 
where w = 0.076 and w\ = 0.527 are estimated. Be- 
yond p c the width decreases as: w{p,L) = \p — p c \~ v 
with v ss 0.82 is estimated. 

The time-displaced autocorrelation of the fluctuating 
mass per site is defined as: 



C(t,p,L) 



(n(t + t,p, L)n(t ,p, L)) - (n(p, L)f 
(n 2 (p,L))-(n(p,L)) 2 



(4) 



This autocorrelation is observed to decay exponentially 
as: C(t,p,L) ~ exp(—t/r(p,L)) where r(p,L) is the cor- 
relation time. On a semi-log plot of C(t,p,L) vs. t the 
slope of the plot gives the value of the correlation time 
t(j>, L) which is measured for different probabilities p and 
for different L values. For a given system size the corre- 
lation time is maximum at p c and then decreases mono- 
tonically with increasing \p — p c \. Also r(p, L) increases 
with L at a given p. A scaling plot of the data collapses 
very nicely as (Fig. 3): 



r(p,L)L 



-3.46 



^((P-Pc)^ 1 - 5 ) 



(5) 



At Pc, t(j? c ,L) increases as L M where /i is estimated to 
be 3.45 ±0.10. 

Fourier transform of the autocorrelation function 
C(t,p, L) is known as the spectral density or power spec- 
trum S(f,p,L) defined as 



S(f,P,L) 



' ift C(t,p,L)dt 



(6) 



In Fig. 4(a) we show the plot of the autocorrelation func- 
tion C(t,p c ,L) with scaled time t/L 2 for a system size 
L — 64 and exactly at p c . A straight line plot on a semi- 
log scale implies an exponential decay of the correlation 
function. In Fig. 4(b) we show the Fourier transform 
of the autocorrelation function plotted in Fig. 4(a) gen- 
erated by the plotting routine i xmgrace\ On a double 
logarithmic scale the power spectrum S(f,p c ,L) vs. the 
frequency / plot gives a very good straight line for the 
intermediate range of frequencies implying a power law 
decay of the spectral density: S(f,p c ,L) ~ f^ 13 . We es- 
timate j3 w 1 showing the existence of 1// type of noise 
in the power spectrum. 

The avalanche size distributions for both particle as 
well as hole avalanches are measured. It is observed that 
in the range of p > p c the particle avalanche sizes are of 
widely varying magnitudes and of all length scales where 
as the hole avalanche sizes are very small and of the order 
of unity. Opposite is the situation for the range p < p c - 
At p c however both the particle as well as hole avalanche 
size distributions are similar and they are found to follow 
a stretched exponential distribution like: 



P(s) ~ exp(— as 7 ) 



(7) 



where 7 is estimated to be around 0.4. Away from this 
critical point p c , the particle avalanches have power law 
distribution for p > p c and hole avalanches follow power 
law distributions for p < p c . Particle avalanche size dis- 
tributions are calculated at p = 0.51 and for system sizes 
L = 256, 512, 1024 and 2048. These distributions are very 
similar to the avalanche distributions in BTW model. For 
small avalanche sizes they do follow a power law distri- 
bution P(s) ~ s~ T where the exponent r slowly varies 
with the system size and gradually increases towards 1.2. 
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The large avalanches have multi-fractal distribution and 
simple scaling does not work for the full distribution 
[ill . Also the average avalanche size, area and life 
times have system size dependances very similar to those 
in the BTW sandpile: (s(L)> ~ L 2 , (a(L)) ~ L 172 and 
(t(L))~L. 

To summarize, we have studied a new sandpile model 
where both particle as well as hole avalanches are created. 
Their relative stengths are tuned by a parameter p vary- 
ing between and 1 which is the probability for adding 
a particle and consequently 1—p being the hole addition 
probability. Specifically at p = 1 the system is identical 
to the ordinary BTW model for only particle avalanches. 
Similarly at p = there is only hole avalanches and their 
distribution are very similar to the avalanche size distri- 
bution for the BTW model. In the range 1/2 < p < 1 
there are particle as well as hole avalanches, but the net 
current is due to the particles which flows into the bulk of 
the system. Critical behavior of the particle avalanches 
are observed to have multi-scaling behavior and is very 
similar to those of the BTW model. Opposite situation 
happens in the range < p < 1/2 where net current is 
due to holes which flows into the bulk of the system. The 
hole avalanche sizes also have multi-scaling distributions 
very similar to the BTW model. 
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